---
title: "Untitled"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(knitr)
library(rlang)
library(readr)
library(tidyverse)
library(stargazer)
library(ggplot2)
library(RCurl)
library(miceadds)
library(plm)
library(gridExtra)
library(lmtest)
library(multiwayvcov)
library(clusterSEs)
library(clubSandwich)
library(sjmisc)
library(ggplot2)
library(gridExtra)
library(grid)
library(ggeffects)
library(estimatr)
library(dotwhisker)
library(lfe)
library(fect)
library(panelView)
library(ggplotify)
library(cobalt)
```

```{r}
#Download data
data <- read.csv("Mass_Shootings_Gun_Policy_Markarian2023.csv")
```

```{r}
ggplot(data , aes(y=lawtotal, x=as.numeric(year), group=state_abb)) + 
  geom_line( size = 1, color = "grey70") +
  geom_smooth(aes(y=lawtotal, x=year, group=control, linetype = control), se = FALSE, size = 2, color = "black")+
  theme_bw()  +
  scale_linetype_manual(values = c( "dashed", "dotted","twodash"), labels = c("Democrats", "Divided", "Republican"), name = "State Control") +
  scale_y_continuous(name="Restrictive Gun Laws")+
  scale_x_continuous(name="Year", limits=c(1991, 2020), breaks = c(seq(1990,2020, by=2))) +
  geom_text(
    aes(label=state_abb),
    data = data %>% filter(year==2020),
    color = "black",
    hjust = 0,
    size = 3
  ) +
  ggtitle("")

```

```{r}
d1 <- select(data, state_full, white_fatalities.tvp.b, year, control)
d2 <- select(data, state_full, nonwhite_fatalities.tvp.b, year, control)
colnames(d1) <- c("state", "fatalities", "year", "control")
colnames(d2) <- c("state", "fatalities", "year", "control")

d1$Race <- "White"
d2$Race<- "Racial and Ethnic Minority"

d3 <- rbind(d1,d2)

ggplot(subset(d3, year>=1991), aes(x=year, y=fatalities, fill=Race)) + 
  geom_bar(stat='identity') + 
  coord_flip() +
  scale_x_continuous(breaks = seq(1990,2020, 5)) +
  labs(x="Year", y= "Mass Shooting Fatalities", title = "") +
  theme_bw() + 
  scale_fill_grey(start = .2, end = .8)
```

```{r}
data$year <- as.factor(data$year)
pdata <- pdata.frame(data, index = c("state_abb", "year"))
```

########################################################################## 
########################### MAIN TEXT ANALYSIS ########################### 
########################################################################## 

```{r}
m1 <- felm(lawtotal.change.prepost ~  white_fatalities.tvp.b + nonwhite_fatalities.tvp.b | state_abb + year | 0 | state_abb, data=data)

m2 <- felm(lawtotal.change.prepost ~  white_fatalities.tvp.b + nonwhite_fatalities.tvp.b  + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int | state_abb + year | 0 | state_abb, data=data)

m3 <- felm(lawtotal.change.prepost ~  white_fatalities.tvp.b + nonwhite_fatalities.tvp.b   + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int + state_abb*as.numeric(year) | state_abb + year | 0 | state_abb, data=data)
```

```{r}
stargazer(m1, m2, m3, covariate.labels = c("White Fatalities", "Racial and Ethnic Minority Fatalities",  "Divided Government" , "Republican Government", "Lagged Total Restrictive Firearm Laws", "Neighbor White Fatalities", "Neighbor Nonwhite Fatalities", "Regular Session" , "Biennium", "Population Black", "Population Latino", "Population Other Race", "Population Median Income"),   dep.var.labels = c("Number of New Restrictive Firearm Laws") ,  omit.labels = c("State Fixed Effects", "Year Fixed Effects", "State Linear Time Trends"), title = "Table XX: Mass Shootings' Effect on State Firearm Laws", type = "html", omit = c("state" , "year", "state*as.numeric(year)" ), star.char = c("+", "*", "**", "***"),  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),  notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"),  notes.append = F, out="table1_aug2023.doc", omit.stat = c("f", "ser"), font.size = "scriptsize")
```

```{r}
m1df <- broom::tidy(m1) %>% 
  filter(term %in% c("white_fatalities.tvp.b" , "nonwhite_fatalities.tvp.b")) %>% 
  mutate(model = "Model 1")
  
m2df <- broom::tidy(m2) %>% 
  filter(term %in% c("white_fatalities.tvp.b" , "nonwhite_fatalities.tvp.b")) %>% 
  mutate(model = "Model 2")

m3df <- broom::tidy(m3) %>% 
  filter(term %in% c("white_fatalities.tvp.b" , "nonwhite_fatalities.tvp.b")) %>% 
  mutate(model = "Model 3")

models <- rbind(m1df, m2df, m3df) %>%                                     # make model variable
    relabel_predictors(c(
        white_fatalities.tvp.b = "Number of\nWhite Fatalities",
        nonwhite_fatalities.tvp.b = "Number of\nREM Fatalities"
    ))

dwplot(models, ci = 0.95,
       vline = geom_vline(
           xintercept = 0,
           colour = "grey60",
           linetype = 2
       ), dot_args = list(aes(shape = model))) +
    theme_bw(base_size = 11) + xlab("Coefficient Estimate") + ylab("") +
    labs(subtitle = "") +
    theme(
        plot.title = element_text(face = "bold"),
        legend.position = "bottom",
        legend.justification = c(0, 0),
        legend.background = element_rect(colour = "grey80"),
        legend.title.align = .5
    ) +
    scale_colour_grey(
        start = .1,
        end = .7,
        breaks = c("Model 1", "Model 2", "Model 3"),
        name = "Model",
        labels = c("TWFE", "TWFE & Controls", "TWFE, Controls, & Linear Trends")
    ) +
    scale_shape_discrete(
          breaks = c("Model 1", "Model 2", "Model 3"),
        name = "Model",
        labels = c("TWFE", "TWFE & Controls", "TWFE, Controls, & Linear Trends")
    ) +
    guides(
        shape = guide_legend(""), 
        colour = guide_legend("")) +

        theme(legend.direction="horizontal")
```

# HHE Robust Estimates
```{r}
#Majority White Mass Shooting
m1 <- fect(lawtotal.change.prepost ~    majority_white_shootings, data = data, index = c("state_abb", "year"), 
  method = "fe", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, placeboTest = TRUE, placebo.period = c(-3, 0), vartype = "jackknife")

p1 <- plot(m1, main = "", ylab = "Effect on Change in Gun Laws", xlab = "Time Since Majority White Mass Shooting", stats = "equiv.p") +
  theme(axis.title = element_text(size = 8))
p2 <- plot(m1, type = "equiv", main = "", ylab = "Effect on Change in Gun Laws", xlab = "Time Since Majority White Mass Shooting", show.points = FALSE)
p3 <- plot(m1, type = "exit", main = "", ylab = "Effect on Change in Gun Laws", xlab = "Time Since Majority White Mass Shooting", show.points = FALSE)

#Majority REM Mass Shooting
m2 <- fect(lawtotal.change.prepost ~    majority_nonwhite_shootings, data = data, index = c("state_abb", "year"), 
  method = "fe", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, placeboTest = TRUE, placebo.period = c(-3, 0), vartype = "jackknife")

p4 <- plot(m2, main = "", ylab = "Effect on Change in Gun Laws", xlab = "Time Since Majority REM Mass Shooting",  proportion =  .45, stats = "equiv.p")
p5 <- plot(m2, type = "equiv", main = "", ylab = "Effect on Change in Gun Laws", xlab = "Time Since Majority REM Mass Shooting", show.points = FALSE)
p6 <- plot(m2, type = "exit", main = "", ylab = "Effect on Change in Gun Laws", xlab = "Time Since Majority White Mass Shooting", show.points = FALSE)

#Graph Edits
p1 +
  theme(axis.title = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12))

p4 +
  theme(axis.title = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.text.x = element_text(size = 12))

#Adding linear time trends
fect(lawtotal.change.prepost ~    majority_white_shootings, data = data, index = c("state_abb", "year"), 
  method = "polynomial", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, degree=1, placeboTest = TRUE, placebo.period = c(-3, 0), vartype = "jackknife")
fect(lawtotal.change.prepost ~    majority_nonwhite_shootings, data = data, index = c("state_abb", "year"), 
  method = "polynomial", force = "two-way", se = TRUE, parallel = TRUE, nboots = 200, degree=1, placeboTest = TRUE, placebo.period = c(-3, 0), vartype = "jackknife")

p2

p5

```

#Demographic Control
```{r}
m1 <- felm(lawtotal.change.prepost ~  white_fatalities.tvp.b + nonwhite_fatalities.tvp.b  + victim_18under + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int | state_abb + year | 0 | state_abb, data=data)

m2 <- felm(lawtotal.change.prepost ~  white_fatalities.tvp.b + nonwhite_fatalities.tvp.b  + female_victim + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int | state_abb + year | 0 | state_abb, data=data)

m3 <- felm(lawtotal.change.prepost ~  white_fatalities.tvp.b + nonwhite_fatalities.tvp.b  + school + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int | state_abb + year | 0 | state_abb, data=data)

m4 <- felm(lawtotal.change.prepost ~  white_fatalities.tvp.b + nonwhite_fatalities.tvp.b + top_tertile_income_FAT  + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int| state_abb + year | 0 | state_abb, data=data)
```

```{r}
stargazer(m1, m2, m3, m4, covariate.labels = c("White Fatalities", "Nonwhite Fatalities",  "Fatalities Age 18 and Under", "Female Fatalities", "School Shooting Fatalities", "Fatalities in Higher Income Census Tracts", "Divided Government" , "Republican Government", "Lagged Total Restrictive Firearm Laws", "Neighbor White Fatalities", "Neighbor Nonwhite Fatalities", "Regular Session" , "Biennium", "Population Black", "Population Latino", "Population Other Race", "Population Median Income"),   dep.var.labels = c("Change in Firearm Laws") ,  omit.labels = c("State Fixed Effects", "Year Fixed Effects", "State Linear Time Trends"), title = "Table XX: Mass Shootings' Effect on State Firearm Laws", type = "html", omit = c("state" , "year", "state*as.numeric(year)" ), star.char = c("+", "*", "**", "***"),  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),  notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"),  notes.append = F, out="tableSIage_aug2023.doc", omit.stat = c("f", "ser"), font.size = "scriptsize")
```

```{r}

  
m1df <- broom::tidy(m1) %>% 
  filter(term %in% c("white_fatalities.tvp.b" , "nonwhite_fatalities.tvp.b", "victim_18under")) %>% 
  mutate(model = "Model 1")

m2df <- broom::tidy(m2) %>% 
  filter(term %in% c("white_fatalities.tvp.b" , "nonwhite_fatalities.tvp.b", "female_victim")) %>% 
  mutate(model = "Model 2")

m3df <- broom::tidy(m3) %>% 
  filter(term %in% c("white_fatalities.tvp.b" , "nonwhite_fatalities.tvp.b", "school")) %>% 
  mutate(model = "Model 3")

m4df <- broom::tidy(m4) %>% 
  filter(term %in% c("white_fatalities.tvp.b" , "nonwhite_fatalities.tvp.b", "top_tertile_income_FAT")) %>% 
  mutate(model = "Model 4")

models <- rbind(m1df, m2df, m3df, m4df) %>%                                     # make model variable
    relabel_predictors(c(
        white_fatalities.tvp.b = "Number of\nWhite Fatalities",
        nonwhite_fatalities.tvp.b = "Number of\nREM Fatalities",
        school = "Number of Fatalities\nin School",
        victim_18under = "Number of Fatalities\nUnder 18",
        female_victim = "Number of Fatalities\nFemale",
        top_tertile_income_FAT = "Number of Fatalities\nin Higher Income\nCensus Tracts"
    ))

g1 <- dwplot(models, ci = 0.95,
       vline = geom_vline(
           xintercept = 0,
           colour = "grey60",
           linetype = 2
       ), dot_args = list(aes(shape = model))) +
    theme_bw(base_size = 11) + xlab("Coefficient Estimate") + ylab("") +
    labs(subtitle = "") +
    theme(
        plot.title = element_text(face = "bold"),
        legend.position = c(0.007, 0.01),
        legend.justification = c(0, 0),
        legend.background = element_rect(colour = "grey80"),
        legend.title.align = .5
    ) +
    scale_colour_grey(
        start = .1,
        end = .7,
        breaks = c("Model 1", "Model 2", "Model 3", "Model 4"),
        name = "Model",
        labels = c("Age Control", "Gender Control", "School Control", "Neigh. Income Control")
    ) +
    scale_shape_discrete(
          breaks = c("Model 1", "Model 2", "Model 3", "Model 4"),
        name = "Model",
        labels = c("Age Control", "Gender Control", "School Control", "Neigh. Income Control")
    ) +
    guides(
        shape = guide_legend(""), 
        colour = guide_legend("")) +
        theme(legend.direction="horizontal",
              legend.position="bottom")


g1
```

##################################################################################
############################ Supplemental Information ############################
##################################################################################

## Race Disaggregated
```{r}
m1 <- felm(lawtotal.change.prepost ~  white_fatalities.tvp.b + black_fatalities.tvp.b + hisp_fatalities.tvp.b + asian_fatalities.tvp.b + other_fatalities.tvp.b + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int | state_abb + year | 0 | state_abb, data=data)

m2 <- felm(lawtotal.change.prepost ~  white_fatalities.tvp.b + black_fatalities.tvp.b + hisp_fatalities.tvp.b + asian_fatalities.tvp.b + other_fatalities.tvp.b + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int + state_abb*as.numeric(year)| state_abb + year | 0 | state_abb, data=data)
```

```{r}
stargazer(m1, m2, covariate.labels = c("White Fatalities", "Black Fatalities", "Hispanic Fatalities", "Asian Fatalities", "Other Fatalities", "Divided Government" , "Republican Government", "Total Firearm Laws", "Neighbor White Fatalities", "Neighbor Nonwhite Fatalities", "Regular Session" , "Biennium", "Population Black", "Population Latino", "Population Other Race", "Population Median Income"),  dep.var.labels = c("Change in State Firearm Laws") , type = "html",   omit.labels = c("State Fixed Effects", "Year Fixed Effects", "State Linear Time Trends"), title = "Table XX: Mass Shootings' Effect on Firearm Laws", omit = c("state" , "year", "state*as.numeric(year)" ), star.char = c("+", "*", "**", "***"),  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),  notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"),  notes.append = F, out="table_racedisagg_aug2023.doc", omit.stat = c("f", "ser"), font.size = "scriptsize")
```

## Proportions test 
```{r }
#Using Lm() function because felm() does not support predict() 
m1 <- lm(lawtotal.change.prepost ~  total_fatalities.tvp*prop.white  + state_abb + as.factor(year), data=data)
se1 <- coeftest(m1, vcov = vcovCL(m1, cluster = data$state_abb, type = "HC1"))

m2 <- lm(lawtotal.change.prepost ~  total_fatalities.tvp*prop.white  + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int + state_abb + year, data=data)
se2 <- coeftest(m2, vcov = vcovCL(m2, cluster = data$state_abb, type = "HC1"))

m3 <- lm(lawtotal.change.prepost ~  total_fatalities.tvp*prop.white  + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int + state_abb + year + state_abb*as.numeric(year), data=data)
se3 <- coeftest(m3, vcov = vcovCL(m3, cluster = data$state_abb, type = "HC1"))
```

```{r}
stargazer(m1, m2, m3,  se = list(se1[,2], se2[,2], se3[,2]),  covariate.labels = c("Fatalities", "Proportion White",  "Divided Government" , "Republican Government", "Total Firearm Laws", "Neighbor White Fatalities", "Neighbor Nonwhite Fatalities", "Regular Session" , "Biennium", "Population Black", "Population Latino", "Population Other Race", "Population Median Income", "Fatalities*Proportion White"),  dep.var.labels = c("Change in State Firearm Laws") , type = "html",   omit.labels = c("State Fixed Effects", "Year Fixed Effects", "State Linear Time Trends"), title = "Table XX: Mass Shootings' Effect on Firearm Laws", omit = c("state" , "year", "state*as.numeric(year)" ), star.char = c("+", "*", "**", "***"),  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),  notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"),  notes.append = F, out="table_prop_aug2023.doc", omit.stat = c("f", "ser"), font.size = "scriptsize")
```

```{r}
df1 <- ggpredict(m1, terms = c("total_fatalities.tvp" , "prop.white [0,1]"), vcov_fun = "vcovCR", vcov_type = "CR0",
  vcov_args = list(cluster = data$state_abb))


ggplot(df1, aes(x, predicted)) + 
    geom_line(aes(color=group)) +
    geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=group, color=group), alpha=0.15, linetype="dashed") + 
  labs(subtitle = "", x = "Fatalities", y= "Change in Gun Laws") + 
  scale_color_manual(values = c( "black", "grey50"), labels = c("Racial and Ethnic\nMinority Fatalities", "White Fatalities")) + 
  scale_fill_manual(values = c("black", "grey50"), labels = c("Racial and Ethnic\nMinority Fatalities", "White Fatalities")) +
     theme_bw()  +
     geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  guides(color=guide_legend(""), fill=guide_legend("")) +
  scale_x_continuous(limits = c(0,30))+
  scale_y_continuous(limits = c(-5,10))
```

## RAND Data
```{r}
m1 <- felm(randrest.change.prepost ~  white_fatalities.tvp.b + nonwhite_fatalities.tvp.b | state_abb + year | 0 | state_abb, data=data)

m2 <- felm(randrest.change.prepost ~  white_fatalities.tvp.b + nonwhite_fatalities.tvp.b  + control + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int | state_abb + year | 0 | state_abb, data=data)

m3 <- felm(randrest.change.prepost ~  white_fatalities.tvp.b + nonwhite_fatalities.tvp.b + control + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int + state_abb*as.numeric(year)| state_abb + year | 0 | state_abb, data=data)
```

```{r}
stargazer(m1, m2, m3, covariate.labels = c("White Fatalities", "Racial and Ethnic Minority Fatalities",  "Divided Government" , "Republican Government", "Neighbor White Fatalities", "Neighbor Nonwhite Fatalities", "Regular Session" , "Biennium", "Population Black", "Population Latino", "Population Other Race", "Population Median Income"),   dep.var.labels = c("Change in Firearm Laws") ,  omit.labels = c("State Fixed Effects", "Year Fixed Effects", "State Linear Time Trends"), title = "Table XX: Mass Shootings' Effect on State Firearm Laws", type = "html", omit = c("state" , "year", "state*as.numeric(year)" ), star.char = c("+", "*", "**", "***"),  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),  notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"),  notes.append = F, out="table_rand_aug2023.doc", omit.stat = c("f", "ser"), font.size = "scriptsize")
```

## PLM Random Effects ~ Both
```{r}
m1 <- plm(lead.lawtotal ~ white_fatalities.tvp.b + nonwhite_fatalities.tvp.b + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int, data=data, model = "within", index = c("state_abb", "year"))
se1 <- coef_test(m1, vcov = "CR1")

m2 <- plm(lead.lawtotal ~ white_fatalities.tvp.b + nonwhite_fatalities.tvp.b + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int, data=data, model = "random", index = c("state_abb", "year"))
se2 <- coef_test(m2, vcov = "CR1")
```

```{r}
stargazer(m1,  m2, se = list( se1[,3],  se2[,3]), covariate.labels = c("White Fatalities", "Nonwhite Fatalities", "Lagged Total Restrictive Firearm Laws", "Divided Government" , "Republican Government",  "Neighbor White Fatalities", "Neighbor Nonwhite Fatalities", "Regular Session" , "Biennium", "Population Black", "Population Latino", "Population Other Race", "Population Median Income"),  dep.var.labels = c( "Lead Firearm Laws") , column.labels = c("Within Effects", "Random Effects"), title = "Table SI-2.1: Mass Shootings' Effect on Proposed State Firearm Laws", type = "html", star.char = c("+", "*", "**", "***"),  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),  notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"),  notes.append = F, out="table_random_aug2023.doc", omit.stat = c("f", "ser"), font.size = "scriptsize")
```


## Cumalitive sum
```{r}
m1 <- felm(lead.lawtotal ~  cum_white_fatalities + cum_nonwhite_fatalities | state_abb + year | 0 | state_abb, data=data)

m2 <- felm(lead.lawtotal ~  cum_white_fatalities + cum_nonwhite_fatalities + control  + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int | state_abb + year | 0 | state_abb, data=data)

m3 <- felm(lead.lawtotal ~  cum_white_fatalities + cum_nonwhite_fatalities + control  + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int + state_abb*as.numeric(year) | state_abb + year | 0 | state_abb, data=data)
```

```{r}
stargazer(m1, m2, m3, covariate.labels = c("Cumulative White Fatalities", "Cumulative Racial and Ethnic Minority Fatalities",  "Divided Government" , "Republican Government", "Neighbor White Fatalities", "Neighbor Nonwhite Fatalities", "Regular Session" , "Biennium", "Population Black", "Population Latino", "Population Other Race", "Population Median Income"),   dep.var.labels = c("Total Firearm Laws t+1") ,  omit.labels = c("State Fixed Effects", "Year Fixed Effects", "State Linear Time Trends"), title = "Table XX: Mass Shootings' Effect on State Firearm Laws", type = "html", omit = c("state" , "year", "state*as.numeric(year)" ), star.char = c("+", "*", "**", "***"),  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),  notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"),  notes.append = F, out="table_cumu_aug2023.doc", omit.stat = c("f", "ser"), font.size = "scriptsize")
```

## Rolling Fatalities Analysis
```{r}
data <- data %>%
  group_by(state_abb) %>%
  arrange(year) %>%
  mutate(rolling_white_fat = white_fatalities.tvp.b + lagged.white_fatalities.tvp.b,
         rolling_rem_fat = nonwhite_fatalities.tvp.b + lagged.nonwhite_fatalities.tvp.b)


m1 <- felm(lawtotal.change ~ rolling_white_fat + rolling_rem_fat | state_abb + year | 0 | state_abb, data=data)

m2 <- felm(lawtotal.change ~  rolling_white_fat + rolling_rem_fat  + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int | state_abb + year | 0 | state_abb, data=data)

m3 <- felm(lawtotal.change ~  rolling_white_fat + rolling_rem_fat   + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int + state_abb*as.numeric(year) | state_abb + year | 0 | state_abb, data=data)
```

```{r}
stargazer(m1, m2, m3, covariate.labels = c("White Fatalities In Last Two Years", "Racial and Ethnic Minority Fatalities In Last Two Years",  "Divided Government" , "Republican Government", "Lagged Total Restrictive Firearm Laws", "Neighbor White Fatalities", "Neighbor Nonwhite Fatalities", "Regular Session" , "Biennium", "Population Black", "Population Latino", "Population Other Race", "Population Median Income"),   dep.var.labels = c("Number of New Restrictive Firearm Laws") ,  omit.labels = c("State Fixed Effects", "Year Fixed Effects", "State Linear Time Trends"), title = "Table XX: Mass Shootings' Effect on State Firearm Laws", type = "html", omit = c("state" , "year", "state*as.numeric(year)" ), star.char = c("+", "*", "**", "***"),  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),  notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"),  notes.append = F, out="table_roling_aug2023.doc", omit.stat = c("f", "ser"), font.size = "scriptsize")
```


## Party Interaction Models
```{r}
#Need to pre-define class or predict() throws an error
data$year <- as.factor(data$year)

#Using Lm() function because felm() does not support predict() 
m1 <- lm(lawtotal.change.prepost ~  white_fatalities.tvp.b*control + nonwhite_fatalities.tvp.b*control +  state_abb + as.factor(year), data=data)
se1 <- coeftest(m1, vcov = vcovCL(m1, cluster = data$state_abb, type = "HC1"))

m2 <- lm(lawtotal.change.prepost ~  white_fatalities.tvp.b*control + nonwhite_fatalities.tvp.b*control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int + state_abb + as.factor(year), data=data)
se2 <- coeftest(m2, vcov = vcovCL(m2, cluster = data$state_abb, type = "HC1"))

m3 <- lm(lawtotal.change.prepost ~  white_fatalities.tvp.b*control + nonwhite_fatalities.tvp.b*control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int + state_abb + as.factor(year) + state_abb*as.numeric(year), data=data)
se3 <- coeftest(m3, vcov = vcovCL(m3, cluster = data$state_abb, type = "HC1"))
```

```{r}
stargazer(m1, m2, m3, se = list(se1[,2], se2[,2], c(se3[1:89,2], NA, se3[90:nrow(se3),2])),  order = c(1,4,91,90,93,92,3,2), covariate.labels = c("White Fatalities", "Nonwhite Fatalities", "White Fatalities*Republican Gov", "Nonwhite Fatalities*Republican Gov", "White Fatalities*Divided Gov","Nonwhite Fatalities*Divided Gov", "Divided Government" , "Republican Government"),   dep.var.labels = c("Change in Firearm Laws") ,  omit.labels = c("State Fixed Effects", "Year Fixed Effects", "State Linear Time Trends"), title = "Table XX: Mass Shootings' Effect on State Firearm Laws", type = "html", omit = c("state" , "year", "state*as.numeric(year)" ), star.char = c("+", "*", "**", "***"),  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),  notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"),  notes.append = F, out="table_partinint_aug2023.doc", omit.stat = c("f", "ser"), font.size = "scriptsize")
```

```{r}
df1 <- ggpredict(m1, terms = c("white_fatalities.tvp.b [0:30, by=10]" ,"control"), vcov_fun = "vcovCR", vcov_type = "CR0",
  vcov.args = list(cluster = data$state_abb))

pd <- position_dodge(width = 1.8)
g1 <- ggplot(df1, aes(x, predicted, color=group)) + 
     geom_line(aes(color=group, linetype=group), size=1.5, alpha=0.5) +
 geom_errorbar(aes(ymin=conf.low, ymax=conf.high), alpha=0.5, position = pd, width=1.5) +
   scale_color_manual(labels = c("Democratic-Control", "Divided Government", "Republican-Control"), values = c("grey30", "black", "grey60")) + 
  scale_linetype_manual(labels = c("Democratic-Control", "Divided Government", "Republican-Control"), values = c("solid","dotted","longdash")) +
      theme_bw() +  
      geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
   guides(color=guide_legend(""),
          linetype=guide_legend("")) +
   scale_x_continuous(limits = c(0,32)) +
  scale_y_continuous() +
  theme(panel.grid.minor = element_blank(),
        legend.direction = "horizontal",
        legend.background = element_rect(colour = "grey80")) +
  labs(subtitle = "", x = "White Fatalities", y= "Change in Gun Laws") +
 scale_y_continuous(limits = c(-13,13))

g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

mylegend<-g_legend(g1)

g1 <- g1 + theme(legend.position="none")
g1
```

```{r}
df2 <- ggpredict(m1, terms = c("nonwhite_fatalities.tvp.b [0:30, by=10]" ,"control"), vcov_fun = "vcovCR", vcov_type = "CR0",
  vcov.args = list(cluster = data$state_abb))

pd <- position_dodge(width = 1.5)
g2 <- ggplot(df2, aes(x, predicted, color=group)) + 
     geom_line(aes(color=group, linetype=group), size=1.5, alpha=0.5) +
 geom_errorbar(aes(ymin=conf.low, ymax=conf.high), alpha=0.5, position = pd, width=1.5) +
   scale_color_manual(labels = c("Democratic-Control", "Divided Government", "Republican-Control"), values = c("grey30", "black", "grey60")) + 
  scale_linetype_manual(labels = c("Democratic-Control", "Divided Government", "Republican-Control"), values = c("solid","dotted","longdash")) +
      theme_bw() +  
      geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
   guides(color=guide_legend(""),
          linetype=guide_legend("")) +
   scale_x_continuous(limits = c(0,32)) +
  scale_y_continuous() +
  theme(panel.grid.minor = element_blank(),
        legend.position = "none",
        legend.box = "true") +
  labs(subtitle = "", x = "REM Fatalities", y= "") +
 scale_y_continuous(limits = c(-13,13))

lay = rbind(
            c(2,2,3,3),
            c(2,2,3,3),
            c(2,2,3,3),
            c(2,2,3,3),
            c(2,2,3,3),
            c(1,1,1,1))

p3 <- grid.arrange(mylegend,g1,g2, layout_matrix=lay)
p3
```

#Shooter Race
```{r}
m1 <- felm(lawtotal.change.prepost ~   white_fatalities.tvp.b + nonwhite_fatalities.tvp.b  + white_shooter  + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int | state_abb + year | 0 | state_abb, data=data)


m2 <- felm(lawtotal.change.prepost ~   white_fatalities.tvp.b + nonwhite_fatalities.tvp.b  + nonwhite_shooter  + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int | state_abb + year | 0 | state_abb, data=data)
```

```{r}
stargazer(m1, m2, covariate.labels = c("White Fatalities", "REM Fatalities", "White Shooter Fatalities", "REM Shooter Fatalities", "Divided Government" , "Republican Government", "Lagged Total Restrictive Firearm Laws", "Neighbor White Fatalities", "Neighbor REM Fatalities", "Regular Session" , "Biennium", "Population Black", "Population Latino", "Population Other Race", "Population Median Income"),   dep.var.labels = c("Change in Firearm Laws") ,  omit.labels = c("State Fixed Effects", "Year Fixed Effects", "State Linear Time Trends"), title = "Table XX: Mass Shootings' Effect on State Firearm Laws", type = "html", omit = c("state" , "year", "state*as.numeric(year)" ), star.char = c("+", "*", "**", "***"),  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),  notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"),  notes.append = F, out="table_shooter_aug2023.doc", omit.stat = c("f", "ser"), font.size = "scriptsize")
```

```{r}
m1 <- felm(lawtotal.change.prepost ~   whiteonwhite + remonwhite + remonrem + whiteonrem + control + lag.lawtotal + ngbr_white_fatalities.tvp.b + ngbr_nonwhite_fatalities.tvp.b  + reg.session.both + biennium + dem_black.int + dem_latino.int + dem_raceother.int + dem_medincome.int | state_abb + year | 0 | state_abb, data=data)
```

```{r}
stargazer(m1, covariate.labels = c("White Fatalities X White Shooter", "White Fatalities X REM Shooter", "REM Fatalities X REM Shooter", "REM Fatalities X REM Shooter", "Divided Government" , "Republican Government", "Lagged Total Restrictive Firearm Laws", "Neighbor White Fatalities", "Neighbor REM Fatalities", "Regular Session" , "Biennium", "Population Black", "Population Latino", "Population Other Race", "Population Median Income"),   dep.var.labels = c("Change in Firearm Laws") ,  omit.labels = c("State Fixed Effects", "Year Fixed Effects", "State Linear Time Trends"), title = "Table XX: Mass Shootings' Effect on State Firearm Laws", type = "html", omit = c("state" , "year", "state*as.numeric(year)" ), star.char = c("+", "*", "**", "***"),  star.cutoffs = c(0.1, 0.05, 0.01, 0.001),  notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"),  notes.append = F, out="table_shooterXfata_aug2023.doc", omit.stat = c("f", "ser"), font.size = "scriptsize")
```

